% section 6.4.1 'robustness' / savings
     % this file shows the differences in CE for GHHW and GLTHI 
     % allowing for savings

clear;
load('./results/allstates_main.mat')
load('./results/model_parameters_main.mat')
addpath functions

MU = model_parameters.MU;
CAPH = model_parameters.CAPH;
gamma = model_parameters.gamma;
rho = model_parameters.rho;
r = model_parameters.r;
T = size(MU,3);
age = [25:(25+T-1)]';

% load objects needed

Pspot_N = simulation_output.Pspot_N;
Ppaid_N = simulation_output.Ppaid_N;
Wio_abi= simulation_output.Wio_abi;
Wio_rea= simulation_output.Wio_rea;
Cg_abi = simulation_output.Cg_abi;
Cg_rea = simulation_output.Cg_rea;
iX = simulation_output.iX;
Nstates = size(Pspot_N,1)-1;
Ppaid_spot = simulation_output.Ppaid_spot;
nInd = size(iX,2);
T = 70;
    
% compute results    
% Ed 13
% maximum allowed level of assets (for grid) under each contract

    amax_phi = 60;
    amax_st = 650;
    [C_PS_abi,A_PS_abi,CE_PS_abi,C_PS_PHI_abi,A_PS_PHI_abi,CE_PS_PHI_abi] = analysis_PS(Wio_abi,CAPH,MU,Ppaid_N,Pspot_N,iX,Nstates,rho,r,gamma,T,amax_phi,amax_st);
    
    % Ed 10
    amax_phi = 40;
    amax_st = 500;

    [C_PS_rea,A_PS_rea,CE_PS_rea,C_PS_PHI_rea,A_PS_PHI_rea,CE_PS_PHI_rea] = analysis_PS(Wio_rea,CAPH,MU,Ppaid_N,Pspot_N,iX,Nstates,rho,r,gamma,T,amax_phi,amax_st);
    
 % graph (doesn't go to paper)

      Cons_PHI_abi =Wio_abi-Ppaid_N;
      Cons_PHI_abi_mean = sum(Cons_PHI_abi.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 
      Cons_PHI_PS_abi_mean = sum(C_PS_PHI_abi.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2)/1000; 
      Cons_HHW_abi =Consumption_HHW(Cg_abi,iX);
      Cons_HHW_abi_mean = sum(Cons_HHW_abi.*(iX~=(Nstates+1)),2)./sum(iX~=(Nstates+1),2); 
      
      
      plot(age,Cons_PHI_PS_abi_mean,age,Cons_PHI_abi_mean,age,Cons_HHW_abi_mean)
      h = legend('GLTHI savings','GLTHI','GHHW');
      set(h,'Location','northwest');
      set(h,'fontsize',14);
      
      saveas(gcf,'./results/PHI_PHI_savs_HHW.png');
    
% table
    % main results to compare  
    CE_GHHW_abi = simulation_output.CE_HHW_abi;
    CE_FB_abi = simulation_output.maintable.CE_FB(2);

    CE_GHHW_rea = simulation_output.CE_HHW_rea;
    CE_FB_rea = simulation_output.maintable.CE_FB(1);

  
    table_precsavings = array2table([CE_FB_rea CE_PS_rea CE_GHHW_rea CE_PS_PHI_rea ;...
                     CE_FB_abi CE_PS_abi CE_GHHW_abi CE_PS_PHI_abi ]);

    table_precsavings = renamevars(table_precsavings,["Var1","Var2","Var3","Var4"],...
                        ["CE_FB","CE_ST","CE_GHHW","CE_GLTHI"]);

    table_precsavings.metric2 = (table_precsavings.CE_GHHW-table_precsavings.CE_GLTHI)./(table_precsavings.CE_GHHW);

    writetable(table_precsavings,'./results/table_precsavings.xlsx');                 